### This script will combine the air temperature data from the HOBO and Air Tables of the microclimate database
### Some unwanted sites are still in the Air raw data, and these will be removed

#Define the in, out, and working directories

in.dir = "/home1/99/jc152199/MicroclimateStatisticalDownscale/MicroclimateData/"
out.dir = "/home1/99/jc152199/MicroclimateStatisticalDownscale/MicroclimateData/"
setwd(in.dir)

# Read in the two datafiles

HOBO = read.csv("DailyHOBOMinTemp.csv", header=T)
iButton = read.csv("DailyAirMinTemp.csv", header=T)

# Generate a unique list of sites from HOBO (this dataset contains only sites that should be included)

sites = unique(HOBO$site)

#Subset data using a for loop and the unique list of sites


t1=NULL
t2=NULL

for (sitex in sites) {

  t1 = iButton[which(iButton$site==sitex),]
  t2 = rbind(t1,t2)
  rm(t1)

  }
  
names(iButton)[5]='micro_min'
names(HOBO)[5]='micro_min'

# Combine two data frames to make one

micro.data = rbind(t2,HOBO)

#Now read in BOM data

in.dir = "/home1/99/jc152199/MicroclimateStatisticalDownscale/BOMData/"
setwd(in.dir)
BOM = read.csv("BOM data daily summary.csv", header =T)

#Create a new data frame with same columns as microclimate data

BOM2 = data.frame(site=BOM$Station_Number, date=as.Date(paste(BOM$Year,BOM$Month,BOM$Day,sep="-"),"%Y-%m-%d"), minairtemp=BOM$MinVal)
BOM2$latdec = as.numeric(BOM2$latdec)
BOM2$longdec = as.numeric(BOM2$longdec)
BOM2$east = as.integer(BOM2$east)
BOM2$north = as.integer(BOM2$north)
BOM2$site = as.factor(BOM2$site)
names(BOM2)[3]='micro_min'

# Change micro.data date field to date format

micro.data$date = as.Date(paste(micro.data$year,micro.data$month,micro.data$day,sep="-"),"%Y-%m-%d")

# Row bind the two data sets

micro = rbind(micro.data[,c(1,6,5)],BOM2)

# 'micro' now contains all the observed microclimate data from CTBCC and BOM sites in the appropriate format
# Now begin working with AWAP data

in.dir = "/home1/99/jc152199/MicroclimateStatisticalDownscale/AWAPData/"
setwd(in.dir)
AWAP = read.csv("AWAP_min_out.csv", header=T)

### AWAP data stored as a single column per site
### Need to remove individual columns (for a single site) along with the date column
### And begin row binding them, the column named after the site must be called 'AWAPmaxtemp' and a site column needs to be added

t.AWAP2=NULL

for (i in 3:40) {

  t.AWAP = data.frame(date=AWAP$date, maxtemp=AWAP[,i], site=names(AWAP)[i])
  t.AWAP2 = rbind(t.AWAP2, t.AWAP)

  }

AWAP = t.AWAP2

# Remove leading 'S' from names of BOM sites in object AWAP

AWAP$site = gsub('S3','3',AWAP$site)

### Dates need to be reformatted to %Y-%m-%d

AWAP$date = as.character(AWAP$date)
AWAP$date= as.Date(paste(substr(AWAP$date,0,4),substr(AWAP$date,5,6),substr(AWAP$date,7,8),sep="-"),"%Y-%m-%d")

# Merge AWAP and Microclimate Datasets

micro.awap2 = merge(micro, AWAP, by=c("site","date"), all.x=T) # This command was used to create original file with max temps

d$date = as.Date(paste(d$date,sep=""),"%Y-%m-%d") # Change date structure of combined micro & AWAP data (original file 'solar_regress_all_days_newtopos_2.csv')

micro.awap2 = merge(d, AWAP, by=c("site","date"), all.x=T) # Merges above file and AWAP min data

names(micro.awap2)[which(names(micro.awap2)=='maxtemp')]='AWAP_min'


#Rename microclimate max temp and AWAP max temp

names(micro.awap2)[which(names(micro.awap2)=='maxtemp')]='AWAP_max'
names(micro.awap2)[which(names(micro.awap2)=='maxairtemp')]='micro_max'

# Change work directory and read in location data for all sites

in.dir = "/home1/99/jc152199/Microclimate Statistical Downscale/Location Data/"
setwd(in.dir)
locs = read.csv("Logger Locs for Downscale.csv", header=T)

# Round location data (UTM to 0 decimals, lats to 4 decimals) and format columns for loc data

locs$east = round(locs$east,0)
locs$north = round(locs$north,0)
locs$latdec = round(locs$latdec,4)
locs$longdec = round(locs$longdec,4)
locs$site = as.factor(locs$site)
locs$east = as.integer(locs$east)
locs$north = as.integer(locs$north)

# Append location data from locs to micro.awap2 (only needs to be done for BOM sites)

sites = unique(micro.awap2$site)

t.micro2=NULL

for (sitex in sites) {

  t.locs = locs[which(locs$site==sitex),]
  
  east = t.locs[1,3]
  
  north = t.locs[1,4]
  
  lat = t.locs[1,5]
  
  long = t.locs[1,6]
  
  t.micro = micro.awap2[which(micro.awap2$site==sitex),]
  
  n = nrow(t.micro)
  
  t.micro$east = rep(east, n)
  
  t.micro$north = rep(north, n)
  
  t.micro$latdec = rep(lat, n)
  
  t.micro$longdec = rep(long, n)
  
  t.micro2 = rbind(t.micro, t.micro2)
    
  }
  
CTBCCmicro_BOM_AWAP = t.micro2

# Remove BOM site 32141, or all rows where AWAP_max is 'NA'

t.micro3 = t.micro2[which(is.na(t.micro2$AWAP_max)==F),]  
   
# Write out a .csv

#write.csv(x=CTBCCmicro_BOM_AWAP, file=paste(out.dir,"Microclimate and AWAP formatted with date.csv",sep=""),row.names=F)
write.csv(x=t.micro3, file=paste(out.dir,"Microclimate and AWAP formatted with date.csv",sep=""),row.names=F)
#write.csv(x=locs, file=paste(out.dir,"Locs.csv", sep=""), row.names=F)

### Fucking around reading in big micro micro max/min topo data

in.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/ToAnalyse/'
setwd(in.dir)

big = read.csv(paste(in.dir,'Micro_AWAPMinMax_Topo.csv',sep=''),header=T)

big$date=as.Date(paste(big$date,sep=""),"%Y-%m-%d")

bigger = merge(big, micro, by=c('site','date'), all.x=F)











